
# This file contains R script that replicates the results reported in 
# Montag, Josef, and James Tremewan. "Let the Punishment Fit the Criminal: An
# Experimental Study." 

# Authors: 

# Josef Montag | josef.montag@gmail.com
# International School of Economics
# Kazakh-British Technical University

# James Tremewan | james.tremewan@wu.ac.at@gmail.com
# Vienna University of Economics and Business, and  
# Vienna Center for Experimental Economics, University of Vienna 

# September 16, 2017 

# The script reproduces Table 1 (summary statistics) and the exact tests
# reported in Section 5.3.

# The main database is "data.Rdata," which is an output of the script
# "Data_Prepare.R" and is located in folder "Data".

# You need to have installed the libraries required by the code.

# The script was tested to work in R version 3.3.3 (2017-03-06), R Core Team.
# 2017. "R: A language and environment for statistical computing." R Foundation
# for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

# Preliminaries {{{1

# Libraries:
library(reshape2)
library(dplyr)
library(tidyr)
library(stargazer)
library(Exact)

# Load data (output of DataPrepare.R):
load('Data/data.Rdata')

# TABLE 1: Summary Statistics {{{1

d.summ <-  d.subject %>%
  select(treatment, donate, randwealth, bonus, wtp, Earnings, Profit, 
         PunOpts1mean, PunOpts2mean, PunOpts3mean, TimeST,
         age, female, crt, Qpolitics) 
summ.stats <- d.summ %>% 
  group_by(treatment) %>%
  summarize_all(funs(mean = mean(., na.rm = T),
                     p25  = fivenum(.)[2],
                     p50  = fivenum(.)[3],
                     p75  = fivenum(.)[4])
  ) %>%
  gather(stat, value, -treatment) %>%
  separate(stat, into = c('var', 'stat'), sep = '_') %>%
  dcast(var ~ treatment + stat)

summ.tests <- d.summ %>% 
  mutate_if(is.logical, as.numeric) %>%
  summarize_each(funs(chisq.t = chisq.test(., treatment)$p.value,
                      mw.t    = wilcox.test(. ~ treatment)$p.value), -treatment
  ) %>%
  gather() %>%
  separate(key, into = c('var', 'stat'), sep = '_') %>%
  dcast(var ~ stat)

tab.summ <- full_join(summ.stats, summ.tests) %>%
  mutate_if(is.numeric, funs(round(., 2))) %>%
  mutate_at(vars(ends_with('mean|.t')), funs(format(., nsmall = 2)))
rownames(tab.summ) <- tab.summ$var
tab.summ <- tab.summ[-1]
tab.summ[grep('PunOpts', rownames(tab.summ)), 
         grep('chisq|mw', colnames(tab.summ))] <- '-'
omit.cls <-  which(colnames(tab.summ) %in% c('chisq.t', 'mw.t'))
tab.summ['TimeST',  -omit.cls] <- 
  paste0(as.integer(tab.summ['TimeST', -omit.cls]) %/% 60,
         ':',
         as.integer(tab.summ['TimeST', -omit.cls]) %% 60)
tab.summ <- rbind(tab.summ, c(sum(-d.subject$Sliders + 1), rep(NA, 3),
                              sum(d.subject$Sliders),      rep(NA, 5)))
tab.names <- c(donate = 'Donation decision',
               bonus = 'Slider Task Bonus', 
               randwealth = 'Random Payment',
               wtp = 'Willingness to Pay',
               Earnings = 'Total Earnings',
               Profit = 'Final Payment',
               PunOpts1mean = 'Avg. punishment: Take All',
               PunOpts2mean = 'Avg. punishment: Take Half',
               PunOpts3mean = 'Avg. punishment: Do Nothing',
               TimeST = 'Time to finish the Slider Task',
               age = 'Age',
               female = 'Female (=1)',
               crt = 'Cognitive Reflection Test (1--3)',
               Qpolitics = 'Political views (0--5)',
               '15' = 'Number of subjects')
tab.summ <- tab.summ[names(tab.names), ]
row.names(tab.summ) <- tab.names
tab.summ[grep('^(Punished|Never|Spiteful)|uniformly', rownames(tab.summ)),
         grep('_p', colnames(tab.summ))] <- NA

# Show the table:
tab.summ

# SECTION 5.3: Exact tests of individual behavior across treatments {{{1

# SUBJECTS WHO NEVER PUNISH: Comparing proportions by donation level {{{2

# Tabulate the proportions:
d.subject %>%
  group_by(donate) %>%
  summarize(mean(NeverPunish)) %>%
  round(2)

# Table for exact tests:
(np <- with(d.subject, table(donate, NeverPunish)))

# Take all vs. take half:
exact.test(np[1:2, ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)
# Take all vs. take nothing:
exact.test(np[c(1, 3), ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)
# Take half vs. take nothing:
exact.test(np[2:3, ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)

# SUBJECTS WHO VARY PUNISHMENT: Comparing proportions by donation level {{{2

# Tabulate the proportions:
d.subject %>%
  group_by(donate) %>%
  summarize(mean(VaryPunishment)) %>%
  round(2)

# Table for exact tests:
(vp <- with(d.subject, table(donate, VaryPunishment)))

# Take all vs. take half:
exact.test(vp[1:2, ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)
# Take all vs. take nothing:
exact.test(vp[c(1, 3), ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)
# Take half vs. take nothing:
exact.test(vp[2:3, ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)

# SPITEFUL PUNISHMENT {{{2

# By treatment

# Tabulate the proportions:
d.subject %>%
  group_by(Sliders) %>%
  summarize(mean(Spite)) %>%
  round(2)

# Table for exact tests:
(sp <- with(d.subject, table(treatment, Spite)))

# Comparing treatments:
exact.test(sp, alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)

# By donation

# Tabulate the proportions:
d.subject %>%
  group_by(donate) %>%
  summarize(mean(Spite)) %>%
  round(2)

# Table for exact tests:
(sp <- with(d.subject, table(donate, Spite)))

# Take all vs. take half:
exact.test(sp[1:2, ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)
# Take all vs. take nothing:
exact.test(sp[c(1, 3), ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)
# Take half vs. take nothing:
exact.test(sp[2:3, ],
           alternative = 'two.sided', method = 'Z-pooled', to.plot = FALSE)

